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Abstract 

We compute the structure factor of the J 1 -J 2 Ising model in an external field on the square lattice within the 
Cluster Variation Method. We use a four point plaquette approximation, which is the minimal one able to 
capture phases with broken orientational order in real space, like the recently reported Ising-nematic phase 
in the model. The analysis of different local maxima in the structure factor allows us to track the different 
phases and phase transitions against temperature and external field. Although the nematic susceptibility 
is not directly related to the structure factor, we show that because of the close relationship between the 
nematic order parameter and the structure factor, the latter shows unambiguous signatures of the presence 
of a nematic phase, in agreement with results from direct minimization of a variational free energy. The 
disorder variety of the model is identified and the possibility that the CVM four point approximation be 
exact on the disorder variety is discussed. 
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1. Introduction 

The structure factor, being a quantity of direct experimental access by neutron scattering and many other 
spectroscopic techniques, is a central quantity in condensed matter physics [1, 2]. Mathematically, it is the 
Fourier transform of the connected pair correlation function, and as such its knowledge gives direct access 
to fluctuations and phase transitions associated to them. Computing the structure factor then amounts 
to compute correlation functions, which is known to be a hard task in statistical mechanics models. In 
order to characterize a phase transition it is often possible to look at simpler one-point quantities, typically 
order parameters, like the magnetization or the density. A qualitative understanding of a phase transition 
can be obtained by simple mean field approximations. If one wants to compute universal quantities, like 
critical exponents, then it is necessary to go beyond mean field approximations, for example through a 
Renormalization Group analysis. But there are special kinds of order which are essentially associated with 
fluctuations and then, even if one is interested in a qualitative description, simple mean field theory does 
not work. This is the case, e.g. of broken orientational phases in systems with competing interactions [3, 
4]. When a competing attraction and repulsion or ferromagnetic and anti-ferromagnetic interactions are 
simultaneously present, the system can develop modulated structures in the form of stripes or bubbles [3]. 
These structures break rotational symmetry of space but may not break translational symmetry, giving rise 
to phases with intermediate (in temperature or external field), purely orientational or nematic-like order, 
in analogy with the nematic phases of liquid crystals [5, 2]. Anisotropic phases with nematic-like order are 
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relevant, e.g. in low dimensional systems like electronic liquid-crystals [ 6 , 7, 8 , 9] and ultrathin ferromagnetic 
films [10, 11, 12, 13]. 

In order to characterize these nematic-like phases from microscopic models it is necessary to go beyond 
naive mean field approximations. In particular, the orientational or nematic order parameter in two dimen¬ 
sional modulated systems is proportional to the difference between correlation functions in two orthogonal 
space directions [4, 14]. As will be better described below, a nematic order parameter can be defined as a 
weighted integral of the structure factor. Then, computing the structure factor gives direct access to the 
nematic-like order in systems with competing interactions and orientational order. 

In a lattice, a systematic way of obtaining better approximations for the thermodynamics of a system is to 
consider clusters of increasing size exactly. The Cluster Variation Method (CVM) is one of a family of cluster 
techniques [15, 16, 17, 18]. Although of mean field character, it allows to improve considerably the locus 
of phase transition lines, specially for systems with competing interactions where naive mean field usually 
gives a very poor approximation to the phase diagram. It is also suitable for computing approximations to 
multipoint correlation functions in a systematic way. The CVM has been applied previously to compute 
the structure factor of a few models as the ferromagnetic Ising model [19] and the two dimensional ANNNI 
model [20]. In reference [21] the authors introduced a general approach for the computation of the structure 
factor within the Cluster Variation Method, and applied it to the Ising model with nearest neighbors (NN), 
next-nearest neighbors (NNN) and plaquette interactions in two and three dimensions. For the case of NN 
and NNN interactions in the square lattice, the so called J 1 -J 2 Ising model, they computed the structure 
factor at zero external field in the paramagnetic phase. The phase transition lines between paramagnetic, 
ferromagnetic and collinear (stripe) phases where characterized and the presence of a disorder line in the 
paramagnetic phase was obtained within the approximation and discussed in relation to the exactly known 
result [22] . Interestingly, in [22] , the four point CVM approximation was proved to render the exact solution 
of the model at zero external field. In a recent work, we applied the CVM to the J 1 -J 2 Ising model in an 
external field [23] and found a nematic phase of the kind discussed above, which had not been identified 
previously. Because of the close relation between the nematic order parameter and the structure factor, we 
decided to extend the method of reference [ 21 ] to compute the structure factor of the model in an external 
field in the whole phase diagram, i.e. also in the relevant ordered phases. 

Results on the J1-J2 Ising model may be relevant to understand part of the phenomenology of high 
temperature superconductors, specially the iron pnictides. For these compounds, a much studied model 
is the quantum Heisenberg J 1 -J 2 [24, 25, 8 ]. This model was shown to have a Ising-nematic phase driven 
by spin fluctuations, which break the Z 4 symmetry of the square lattice, without the development of anti¬ 
ferromagnetic order [24]. Strong spin fluctuations in this 2D system induce a biquadratic or quadrupolar 
interaction leading to Ising-like behavior in spin space and eventually to the presence of an Ising-nematic 
phase. Nevertheless, it is not clear if the quadrupolar coupling is strong enough to apply to the experimental 
compounds which show J1-J2 behavior. Another route to nematic order in the pnictides seems to be 
related with doping. Recent results of Monte Carlo simulations on a model with magnetic, electronic and 
orbital degrees of freedom imply that the nematic phase is enhanced through Fe substitution by impurities, 
i.e. by introducing quenched disorder and magnetic dilution in the parent compound [26]. Very recently, 
Kitada et al. [27] reported on an extensive series of experiments on the layered perovskite RbLaNb 2 07 
transformed by substitution of the Rb on the oxyhalides (MCl)LaNb 207 (M=Mn, Cr, Co), which are two 
dimensional antiferromagnets. Neutron diffraction measurements show the presence of magnetic modulations 
with wave vectors ( 0 , 7 r) for the samples with Mn and Co. The samples with Cr showed instead a ( 7 r, 7 r) 
Neel antiferromagnetic structure. Interestingly, hysteresis measurements on the (CoCl)LaNb 207 compound 
indicate a way to saturation in two steps, as the field is raised. This is interpreted as Ising-like behavior, 
in which the striped ground stated is destabilized by a ferromagnetic component by first flipping half of 
the antiferromagnetic stripes and at a higher field value the other half is flipped, leading to the completely 
saturated state. If confirmed, this is the first compound to show a phenomenology typical of the Ising Ji — J 2 
model studied in the present work. 

In the following, we make a brief discussion of known results on the J1-J2 Ising model, the CVM approach, 
and compute the structure factor of the model in presence of an external field in the whole parameter range. 
We interpret the results in connection with the recently published phase diagram [23] . We also identify the 
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disorder variety of the model in an external field and discuss the possible exactness of the four point CVM 
approximation on this variety. 

2. J 1 -J 2 Ising model and the CVM approximation 

The J 1 -J 2 Ising model on the square lattice is defined by the Hamiltonian: 

n = JiY. s^Sy + J 2 YI ^ (1) 

(xy) {{xy)) ^ 

where {S'a; = ±l,x = l...iV} are N Ising spin variables and h is an external field, {xy) denotes a sum over 
pairs of nearest-neighbors and {{xy)) a sum over pairs of next-nearest-neighbors. In this work we consider 
Ji < 0 and J 2 > 0 representing ferromagnetic NN and anti-ferromagnetic NNN interactions respectively. 
The competition ratio is defined hy n = > 0. 

At zero external field the ground state of the model is ferromagnetic for k < 1/2 and striped (super¬ 
antiferromagnetic) if At > 1/2. The stripe phase has fourfold degeneracy as shown in Figure 1. For h = 0 
the model has been extensively studied [28, 29, 21, 30, 31, 32, 33, 34, 35]. The nature of the thermal phase 
transition from the stripes to a disordered phase for k > 1/2 was controversial. In the most recent studies 
combining Monte Carlo simulations and a series of analytic techniques it has been established that the line 
of phase transitions in the temperature versus k plane is first order for 1/2 < k < 0.67 and is continuous 
with Ashkin-Teller critical behavior for k > 0.67. The critical exponents change continuously in this regime 
between the 4-state Potts model behavior at k = 0.67 to standard Ising criticality for k —>• 00 [33, 34]. 



Figure 1: Sketch of the striped ground state configurations of the J 1 -J 2 model in the square lattice. 

For K > ^ and small magnetic fields the ground state is striped. Stripe order is eventually destroyed 
at a critical field value he = ±2(Ji -|- 2 J 2 ). For h > he all spins are aligned with the field and the ground 
state corresponds to a saturated paramagnet. In a recent work it was found a new equilibrium phase at 
intermediate fields in the h vs T plane, a Ising-nematic phase with uniform magnetization but different 
nearest-neighbor correlations along the two directions of the square lattice, breaking the fourfold rotational 
symmetry [23]. At lower fields a second transition takes place to a full stripe phase with broken rotational 
as well as translational symmetries. In order to admit a nematic phase the system must have enhanced 
fluctuations, this is the reason behind the Ising-nematic phase originally reported in the Heisenberg J 1 -J 2 
model [24]. In the 2D Heisenberg model strong fluctuations are due to the continuous rotational symmetry 
in spin space. In the Ising model this mechanism is absent, but enhanced fluctuations can be induced by 
the switching of a magnetic field. In fact, in a restricted temperature window which induces temperature 
fluctuations, the net effect of an external field on the stripe ground state is to favor one of the Ising directions 
and to weaken the other. For a suitable intensity the field will be responsible for inducing defects on the 
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stripe pattern and eventually an instability leading to the loss of anti-ferromagnetic positional order, but 
still preserving a preferred orientation for the remaining stripe pattern. This is precisely the signature of 
the Ising-nematic phase, a phase with broken positional order (or translation symmetry) but still having 
orientational order of the stripe pattern on the lattice. 

The Cluster Variation Method consists in improving a mean field approximation in a systematic way by 
summing all the degrees of freedom within clusters of size n in an exact way. It amounts to extremize the 
variational free energy: 

Ft = Tr{ptH) + kBTTr{pt\npt), (2) 


where Tr means a trace or a sum over all the relevant degrees of freedom of the Hamiltonian H, and pt is 
a trial density matrix which satisfies the normalization constraint Tr pt = 1. Details on the method can be 
found in the large literature on the subject (see e.g. references [15, 36, 17, 18]). 

In the case of a system with Ising spins {Si = ±1}, the reduced density matrix for a cluster of size n, 
p[^\ can be written as [21]: 


Pi 


(n) 


= 2 “ 


1 + O'kCk 


( 3 ) 


where the sum runs over all sub-clusters with k sites within cluster n , ak = Iliefc k-point 

(k) 

correlation functions are defined hy ({k = Tr akp\ ■ In this way, instead of optimizing with respect to the 
reduced densities, the variational parameters are the k-point correlations Qk which must satisfy: 



( 4 ) 


A hierarchy of approximations to the free energy can be constructed in this way. The simplest one cor¬ 
responds to the 1-point approximation for the density matrices, the usual mean field approximation. The 
2-point approximation is usually called Bethe-Peierls approximation [37, 18]. As discussed in [23] the minimal 
cluster able to capture the emergence of anisotropic nearest-neighbor correlations or spontaneous rotational 
symmetry breaking is the four-point or square approximation. 

Define the correlation functions: 




^yxw 


{Sx) = Tr{Sa;Px) 

{SxSy) = Tr [SxSyP(xy)) 

(SxSz) = Tr {SxS^p^i^xz))) 

(SySxS^) = Tr 

{SySxSw P[yxw]) 

{SxSyS^Sxu) = Tr (SxSyS^Su,pxu^\ , 


( 5 ) 


where the sums over x, {xy), {{xy)), [xyz], refer to all sites, NN pairs, NNN pairs, clusters of three 
sites and squares respectively. They are related to the reduced density matrices by: 


Px 

= 

i (1 + TTIj^Sj:) 

P{xy) 

= 

— (1 -|- TfixSx "^ySy IxySxSy) 

P{{xz)) 

= 

— (1 + TfixSx TfizSz C-xzSxSz) 


= 

^(1 + rUxSx + 'f^ySy + rrizSz + m-ujS-uj 


-b 

^XwSxSyj ~\~ IwzSyjSz ^ZySzSy “h IxySxSy 


-b 

^XzSxSz C-ywSySix kyxwSySxSyj 


-b 

^XWzSxSyjSz ~\~ kixzySyjSzSy “t” ^zyxSzSySx 


-b 

dxyZwSxSyS zSyj^ . 
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( 6 ) 




Then, the free energy of the J 1 -J 2 model in the CVM square approximation can be written as [21, 23]: 


F — J\ ^ ^ Ixy J 2 ^ ^ Cxy h ^ ^ TTlx 

{xy) {{xy)) X 

H“ ^ ^ TT {^Px^OQPx^ xy )logp{ 

xy)) 


X 


{xy) 




^Tr (p»n»;o(7pxn») 


( 7 ) 


After computing the traces one is left with an expression for the variational free energy in terms of a set of 
correlation functions representative of the approximation considered. The form of equation (7) makes clear 
that up to the pair approximation the two directions in the square lattice enter in a completely symmetric 
way, the different pairs of sites are decoupled. It is in the last term of (7), when square plaquettes are 
considered, that the coupling between different directions in space can lead to novel behavior. 

3. The structure factor and the nematic order parameter 

In order to detect orientational order it is natural to define orientational or nematic order parameters. 
The nematic order parameter was introduced originally in the study of ordered phases of liquid crystals [5, 2]. 
More recently, nematic-like order was found to be useful to characterize orientation of interfaces or space 
modulations of some physical density, like electron density in the so-called ” electronic liquid-crystals” [6] or 
spin density in magnetic systems [10, 13]. In these systems a nematic order parameter can be defined, in 
analogy with the nematic order parameter of liquid crystals, as a second-rank symmetric traceless tensor 
which encodes the 180° rotational symmetry of nematic phases. In two dimensions the nematic tensor has 
only one independent entry which can be written as [4]: 



( 8 ) 


where k = (kx, ky) and 6 is the angle between the local wave vector k and a fixed direction in the plane. It 
is clear that the nematic order parameter amounts to compute a weighted average of the structure factor of 
the system. The weighting factor cos ( 29 ) has exactly the symmetry of the nematic phase and then Q will be 
zero if S{k) is isotropic. Then, the nematic order parameter amounts to compute the degree of anisotropy of 
the structure factor of the system. Also, because the structure factor is a quantity of primary experimental 
relevance, it is interesting to be able to compute it in the context of the Cluster Variation Method. This has 
been done in a few previous works [19, 20, 21]. In reference [21] the authors introduced a general method for 
the computation of the structure factor at any level of approximation in the CVM and computed it for the 
paramagnetic phase of the J 1 -J 2 model at zero external field in the four point approximation in the square 
and simple cubic lattices. The starting point is the computation of the two-point connected correlation 
function: 



Then, in the square lattice, the structure factor is simply the discrete Fourier transform of C'c(r): 


L-l 


S{k)= 


( 10 ) 


ri,r2=0 


In this work, we have extended the method introduced in reference [21] in order to compute S{k) also in 
the ordered phases of the J 1 -J 2 model and with the inclusion of a uniform external field h. Details of 
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the calculations are presented in the Appendix. We are interested in particular in identifying signatures 
of the nematic phase in the structure factor. The computation in this case, although straightforward, 
is considerably more cumbersome than in the zero field case for the paramagnetic phase, because of the 
presence of a non zero magnetization and the need to consider the symmetry of the ordered phases. With 
the aim of searching for purely orientational nematic-like phases, i.e. phases without translational order, 
we minimized the CVM free energy of Eq. (7) with the following symmetry in the parameters (related to 
the four points of the elementary square yO'^): rux = my, = m^, Ixw = lyz, Ixy, Iwz, c, kyxw = kzyx, 
kxwz = kwzy and d. This choice implies possible orientational order along the xy or vertical direction. Note 
that local magnetizations on horizontal NN sites are allowed to be different in sign and also in absolute value. 
Correspondingly, the NN correlation functions Irs may be different not only between the horizontal and 
vertical directions but also between the two vertical ones. With these choices the values of NNN correlations 
c and square correlations d are unique. The values of Irs, Crs, krst and d with r,s,t = x,y,z,w that 
minimize the variational free energy (7) were calculated numerically for different temperatures (T) and 
external fields {h) as described in [23]. 

After the symmetry considerations, the structure factor is found to take the form (see the Appendix for 
details): 


S{k) ^ = 2{'^xx+lww) + '&lxwCOSki+‘i{'^xy+lzw)cOSk2 

+ 87a;z [C0s(fci -I- fc 2 ) -I- COs(fci - ^ 2 )] (11) 

where k = (fci,fc 2 ) and the coefficients ■jrs are the inverse pair connected correlations as defined in the 
Appendix. They depend on mr, Irs, Crs- krst and d with r,s,t = x, y, z, w as in equation (A. 8 ). 


4. Results 

4-1. Generalized suseeptibilities and phase transitions 

In this section we describe the evolution of the structure factor with decreasing magnetic field for k = 0.6 
at a fixed reduced temperature T/| Jij = 0.67. Along this line the system goes through two successive phase 
transitions as the external field is lowered: one from the isotropic to the Ising-nematic phase and another 
one, at lower field value, from the nematic to the stripe phase. The structure factor shows a series of 
extremes in the {ki,k 2 ) plane, with heights which evolve with the temperature and magnetic field. These 
extremes are located at wave-vectors (0, 0), (0, ± 7 r) and (± 7 r, 0). The heights of the extremes correspond to 
generalized susceptibilities: 

X{kuk2)=PS{kuk2). (12) 

Typically, at a second order phase transition the generalized susceptibility displays critical behavior at 
particular values of the wave vector, diverging at the critical point. At a first order phase transitions it suffers 
a discontinuous jump. The susceptibilities corresponding to the five extremes are shown in Figure 2. For 
large fields h/\Ji\ > 0.45 the system is in a homogeneous paramagnetic state. In this region jxx = Iww = 7 , 
Ixw = Ixy = Izw = 7 i and Xxz = 72 and the structure factor has the same form as in the zero field case [ 21 ]. 
For h/jJij > 0.47 the maximum corresponds to the peak at the origin k = (0,0) (green line). The other 
extremes (saddles, red and black lines) are located on the axes at the border of the first Brioullin zone, 
k = (± 7 r, 0 ) and k = ( 0 ,± 7 r), have equal heights and cross the green line at hi/] Ji] = 0.47. This point is a 
“Lifshitz point”, where the solutions at non zero wave-vector become metastable. At this point the inverse 
pair correlations satisfy 71 = — 272 . 

A density plot of the structure factor for h/\Ji\ = 0.5 is shown in Fig. 3a. The second remarkable fact is 
that the absolute maxima for h < ior k = (± 7 r, 0) and k = (0, ± 7 r) have equal height until hc 2 /\Ji\ = 0.45 
where they bifurcate. In the whole sector h/\Ji\ > 0.45 the spatial distribution of the magnetization has 
Z4 symmetry, characteristic of the square lattice. At hc 2 /\Ji\ = 0.45 there is a spontaneous breaking 0/Z4 
symmetry, to a phase with a lower, Z 2 symmetry. This phase transition, which was shown to be continuous in 
[23], is a paramagnetic to Ising nematic phase transition. Note that the generalized susceptibilities represent 
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Figure 2: (Color online) Generalized susceptibilities x(fci, ^ 2 ) at fc = (0,0), (0 ± tt) and (±7r, 0) for T /\Ji\ = 
0.67 as function of external field. Inset: zoom around the highest peak in logarithmic scale. 


fluctuations of the magnetization and not of the nematic order parameter and then, although the transition 
is second order, it is not accompanied by a divergence of the magnetic susceptibility. In order to study 
the divergence of the nematic susceptibility it would be necessary to go beyond two point correlations and 
consider their own fluctuations, i.e. four-point correlations. A density plot of the structure factor in the 
Ising-nematic phase for h/\Ji \ = 0.44 is shown in Fig. 3b. 

As can be seen in Fig. 2 the secondary peaks at fc = (0, ±7r) decrease rapidly for h < hc 2 , while the 
peak at the origin grows steadily although at a much lower rate than the peaks at (±7r,0). Upon further 
lowering the field a second singularity appears at a field value hcil\Ji\ = 0.382 (see the inset in Fig. 2). At 
this point the susceptibility at (±7r, 0) changes discontinuously while the peak at (0, ±7r) is negligibly small. 
This transition is not accompanied by a symmetry breaking, instead it is a first order transition in which 
the stripe (translational) order parameter jumps from zero to a finite value for h < hci, as originally found 
and discussed in [23]. In Fig. 3c we show a density plot of the structure factor just below the critical field 
h/\Ji\ = 0.38 < /ici/|Ji|. Our results correspond to the expectations for an Ising-nematic phase, see e.g. 
figure 2 in [ 8 ] and figures 7 and 8 in [26]. 

Disorder line 

It is well known that many one and two dimensional frustrated systems have a disorder variety in 
parameter space which crosses the disordered phase. At one side of the disorder variety the pair correlations 
show a monotonic exponential decay with distance, while on the other side the correlations show a modulation 
or a sinusoidal decay, typical of the presence of frustrated or competing interactions [38, 39, 40, 41]. On the 
disorder variety the pair correlation function usually factorizes along two perpendicular directions C{rx,ry) = 
g{Tx)g{'fy), showing typical one dimensional behavior which in turn allows to obtain, in many cases, the 
exact solution of the model. This is the case of the J 1 -J 2 Ising model on the square lattice. For the zero 
field case, it was shown that the four-point CVM approximation yields the exact solution of the model 
on the disorder variety [22]. In particular, g{u) = exp{—with ^ = —l/ln(Z), I being the nearest- 
neighbor correlation. This implies that the structure factor has a similar one dimensional factorization 
S{kx,ky) = Si{kx)Si{ky), where Si{ki) = sinh(l/^)/[cosh(1/^) — cosfci]. Given this simple behavior, it is 
easy to show that on the disorder variety the inverse pair correlations satisfy the constraint 7^ = —772. 
This is an exact relation for the zero field case which, when written in terms of T/\Ji\ and J2/IJ1I, leads 
to the disorder line shown in Figure 4. When an external field is present the general form of the structure 
factor is given by (11). Nevertheless, one expects that the disorder line should continue to exist in the 
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(c) h/|Ji| = 0.38 (deep in the stripe phase). 

Figure 3: (Color online) Density plot of the structure factor for T/| Ji| = 0.67 and different reduced magnetic 
fields. 


disordered region with k < 0.5, where there is a unique paramagnetic phase with finite magnetization. In 
this context, it is natural to expect that the structure factor should have the same simple form of the h = 0 
case, with only three relevant parameters 7 , 71 and 72 . Then, assuming a one dimensional factorization on 
the disorder variety, as in the h = 0 case, the same relation = —772 defines a disorder surface, where 
in this case each 7 is a function of r/|Ji|, J 2 /\Ji\ and h/\Ji\. The similar form of the structure factor in 
the disordered phases of the zero field and finite field cases suggests that also for finite h the relation which 
defines the disorder variety should be exact. The existence of an exact solution for a system which displays a 
nematic-like phase, even limited to the disordered region of the phase diagram, is an interesting possibility. 
This is an open question which deserves further study and is beyond the scope of the present work. 


8 











Figure 4: (Color online) Phase diagram for h = 0 from the four point approximation of the CVM. The exact 
disorder line is also shown. 


5. Conclusions 

We have computed the pair correlations and structure factor of the J 1 -J 2 square lattice Ising model in 
an external field within the Cluster Variation Method. Our motivation was the analysis of the recently 
reported Ising-nematic phase in a sector of the T — h phase diagram of the model. Because the nematic 
order parameter is a function of the structure factor then the latter should show clear signatures of the 
presence of a nematic-like phase. Considering rather general symmetry conditions for the values of the local 
magnetizations in an elementary square of the lattice, we applied the four point approximation in the Cluster 
Variation Method, which is the minimal approximation capable of detecting the presence of an orientational 
phase of nematic character. We showed that the results for the structure factor are in agreement with the 
corresponding phase diagram reported in [23] and that, although the nematic susceptibility is not directly 
related to the structure factor, the presence of the Ising-nematic phase is clearly evidenced in its behavior. 

It was shown that the disorder variety of the model is defined by a constraint between the inverse 
correlations in an elementary square and that the form of the constraint is the same for the cases with or 
without external field. This, together with the exactness of the four point approximation on the disorder 
variety at zero field makes it plausible that also for finite field the approximation should be exact on the 
disorder variety, a point that deserves further study. 

The observation of nematic-like phases in condensed matter systems, like ultrathin ferromagnetic films 
and electronic liquids, is growing rapidly in recent years and experimental determination of structure factors 
and Fermi surfaces in fermionic systems is being increasingly reported in studies of low dimensional magnetic 
systems at the nanoscale and high temperature superconductor systems, to cite two important examples. 
Then, the analytic determination of the structure factor in suitable approximations like the ones accessible 
within the Cluster Variation Method is a valuable tool to compare with computer simulation studies and 
experimental results on these systems. 


Acknowledgment 

A.G.D. and D.A.S. acknowledge partial financial support by Conselho Nacional de Desenvolvimento 
Cientifico e Tecnologico (CNPq), Brazil fellowship 303917/2013-0. 


9 







Appendix A. Connected correlations within the four-point approximation 


The general method to compute pair correlation functions in the CVM has been studied by some authors 
[19, 20, 42, 21], The method consists in introducing an external magnetic field in the variational potential 
and expressing the correlation functions in terms of successive derivatives of the free energy with respect to 
the magnetic field, as in equation (9). However, it is easier to compute the inverse connected correlation 
function: 



(AT) 


In the four-site approximation, the expression for the inverse correlation function can be obtained from the 
state equation for the magnetization, differentiating the variational free energy (7) with respect to rrii: 


phi 


Y. Tr{SPogp(,y)) + ^ Tr{Sdogp.^a^) 
(iy) 

^Tr{SPogpi), 


and then differentiating with respect to rrij: 


(A.2) 
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dhj 

drrij 


1 Si dpi 1 Si dpij 

2 Pi drrij 4 ^ pij [drrij 


dpij dlij 
dlij dnrij 
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dpi^-u, p)i. p)i 

H- - -^- - - [ 

dliw drrij dlzw drrij 

^ dkzn., ^ dpin^ dhjz 

dkzwi drrij dkijz drrij 


dPinr dljz ^ dc,z ^ dc,^ 

dljz drrij dciz drrij dcjw drrij 

^ ^P'P^ dkjzv. ^ dd^ 

dkjzw dTfij ddijzw dfrij 


dpi.n^ dpi.n^ 

3 ^ _j_ 3 ^ 

drrij dlij 

^Ppf dk^,j 

dkwij drrij 


dhj 

drrij 


(A.3) 


In the general case, all the derivatives different of zero. To compute 

the derivatives is not an easy task. One way, as is mentioned in reference [21], is by differentiating the 
variational free energy with respect to Ixy, Cxz, kxyz and dxyzw- 


dJ^ 

dixy 

OT 

dcxz 


= 0 
= 0 


dJ^ 

dkxyz 


= 0 


dr 

ddxyz 


0 , 


(A.4) 


and then, with respect to rrij, j = x, y, z, w: 

^(d^]=n a ( ^ - 0 

dnij ydlxy J drrij ydkxyz J 

d^ (afx^zu,) ^ 

(A.5) 


Giving the symmetry considerations in the Ji — J 2 model: rrix = rriy, rriw = ruz, Ixw = lyz, Ixy, Iwz, c, 
kyxw = kzyx, kxwz = kwzy and d, we found a system of twenty two linear equations and twenty two variables 
Zi, ..., Zg, Cl,..., C4, fci,..., fcg and dip2, where: 
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h — 
h = 

Cl = 


fci = 
fcs = 
di = 


9lxy 

dirix 

h = 

aix,„ 

drrix 

h = 

dlz.u, 

drrix 

U = 

dlyz 

drrix 

dixy 

drrixtj 

Iq = 

dlxw 

dm-uj 

^7 = 

dlzw 

drriw 

h = 

dlyz 

dmy, 

dCxz 

drrix 

C 2 = 

dCyuj 

drrix 

C 3 = 

dCxz 

dm-w 

C4 = 

dc-uw 

drriy, 

dkyzw 

drrix 

k2 = 

dkxyz 

drrix 

fcs = 

dkz-wx 

drrix 

ki = 

dk-ujxy 

drrix 

dkyzw 

dm-w 

ke = 

dkxyz 

dm-w 

k^ = 

dkz-wx 

dm-u) 

ks = 

dkyjxy 

dmy, 

ddxyzw 

drrix 

d2 = 

ddxyzw 

dm-w 




(A. 


Differentiating (A.2) with respect to rrij, with i,j nearest neighbors, we found that derivatives of the 
kind and where j ^ appear. Following the procedure described before, we 

found a set of twenty two linear equations and twenty two variables l'^, Zg, c'^,..., C 4 , fcg and 

similar to (A. 6 ) with rrij external to the plaquette. All these derivatives are zero. Collecting all the pieces, 
from (A.3), the inverse correlation functions are given by: 




dhj 

dm-. 


= 


^xx A ^yy “t“ '^zz “f ^ww 
^'yxw 

‘^'Ixy “t“ 2y zw 
'^Ixz 

0 


I =J 

{ij)h 

{ij)v 

m) 

otherwise, 


(A.7) 


where the derivatives /3 are the inverse self correlations (i = j), NN inverse correlations in the 

horizontal and vertical {{ij),u) direction, and inverse correlation for NNN {{{ij))), respectively. All 

other correlation functions are zero. 

The coefficients 7 ^ are given by: 
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with: 


drrix 

T7 -Ri ~ ■577(^5 + Qe + Qih + Q212) + TTj {Pio + Pih + P2 h + Ps h 

lo 62 lo^ 

+-P7 14 + P2 Cl + -P7 C2 + Pg k\ + -P5 A:2 + Pi + Pq k^ + -Pg c^i) 

^ dhyj 
dlTlw 

TT R2 — ■ 57 r(Q 6 + Q7 + Q3 ^6 + Qi h) + (-P10 + P7 ^5 + Pi Ig + P2 h 

16 32 16^ 


+P8 Is + Ps C3 + Pi C4 + Pg kg + Pg kg + Pg kj + P3 kg + P7 dg) 


= /3 


dhx 

druu 


Ixy — 


~ 2 2 + -Pi ^5 + P2 Ig + Psh + P7 Is + P2 C3 

+P7 C4 Pg kg + Pg kg + P4 ^7 + Pg kg + Ps ^2) 

^ dhoo 
dmy 

~ 2 2 + Pi + P7 ^2 + Ps ^3 + P2 ^4 + P7 Cl 

+P2 C2 + P4 /ci + Pg kg + Pg fcg + Pg ki + Pg di) 


= /3 


9/lx 

dtriz 

1 1 


— 2 + Pi /s + P7 /e + Ps /7 + P2 /s + P2 C4 

+P7 C3 + P4 /cg + Pg fcg + Pg ^7 + Pg kg + Pg dg) , 


Pi = 


Pg = 


tr 

tr 

tr 

tr 


i^pxnj, 

pxnj, 

Sn, 

pxn» 

s, 

pxn» 


P 3 = 
P 4 = 
Pg = 
Pg = 


tr 


I pxnj, 


= tr 


tr 


tr 


Sn,S^ 

pxfj-w 


SxS^ 

pxnju 


tr 


SxSr. 

^pxn» 

I pho™ 


P7 = 


Pg = 


tr 


= tr 


tr 


= tr 


Pg = tr 


SxSySx 

pxnu, 


Sy,SxSy 

px\^xu 


SySxSu 

pxntp 


S^Sy.Sx 

pxntp 


SxSySxSy 

pX\^‘W 


Pio = 


tr 


px\J-w 


(A.8) 
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Finally, the discrete Fourier transform of the inverse correlation is; 


ri,r2—0 ^ 

= ^{ixx + Iww) + ^Ixw cos ki + + Jzw) cos ^2 

+ 8(7xz) [cos (fci + ^2) + cos (fci - fc2)] . 


(A.9) 
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